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SUMMARY 


A computer code named DOWN has been created to implement a "flat wake theory" for 
the calculation of rotor inflow and wake velocities. The theory was developed by V.E. Baskin of 
the USSR. The code was developed at Princeton University under a Grant from the National 
Aeronautics and Space Administration (NASA). A brief description of the code methodology 
and instructions for its use are given. The code will be available from NASA’s Computer 
Software Management and Information Center (COSMIC). 

PROBLEM DEFINITION 

The prediction of inflow to a helicopter rotor and the wake velocities below and behind 
are vital to the calculation of airloads on the rotor blades of a helicopter and the other com- 
ponents. Inflow and wake velocities affect the local angles-of-attack, and therefore the airloads, 
that the helicopter rotor and other components experience. 

Many theories and analyses have been developed for prediction of both mean and time- 
varying inflow/wake patterns. Computer codes to implement these predictions have grown in 
complexity as the theories have grown in complexity. Rotor performance codes use inflow 
models that range from the assumption of uniform inflow to complex, time-varying, vortex 
filament, "free-wake" models. Even though computer processing time and associated costs have 
dramatically decreased in recent years, there is still need for codes which are simple, fast, and 
accurate for preliminary design and analysis. 

PROBLEM SOLUTION 

A relatively simple analysis for predicting helicopter inf low/ wake velocities called "flat 
wake theory" was implemented in a computer code, named DOWN. The analysis essentially 
treats the rotor wake geometry as rigid without interaction between induced velocities and wake 
structure. The flat wake theory for rotors was developed by V.E. Baskin of the USSR and is 
described in reference 1. An analytic scheme for implementing the theory is described in 
reference 2. The theory and associated assumptions and limitations are described in references 1, 
2 and 3. 

The computer coding of the scheme was implemented by Professors Howard C. Curtiss 
and Robert M. McKillip of Princeton University under National Aeronautics and Space Ad- 
ministration Research Grant No. NAG- 1-1038 (Studies of Rotor Inflow Using a Flat Wake 
Theory Including Correlation with Experiment). The code will predict the three orthogonal in- 
cremental components of flow velocity (as ratios of incremental velocities to tip speed) at any 
point in any plane parallel or perpendicular to the rotor disk. The input to the code is quite 
simple and is entered interactively through the computer keyboard. 

The predictive capability of the coded version of the theory has been correlated with 
flow velocity data of several sources. In general, the coded version of flat wake theory provides 
vertical inflow patterns similar to experimental patterns for helicopter flight speeds greater than 



approximately 60 knots (or, as indicated in reference 1, for rotor advance ratio greater than 
0.15). 


SYMBOLS 


Symbols used in the code are bracketed. 


C T 

(CT) 

rotor thrust coefficient, T /[p a (OR) 2 7rR 2 ] 

G 


circulation "normalizing factor" 

r 


rotor blade station, ft 

R 


rotor radius, ft 

T 


rotor thrust, lbs 

v o 


momentum theory value of average induced velocity, fps 

v o 


v Q normalized with tip speed i.e. v Q /V T 

V 


wind or flight speed, ft/sec 

V t 


rotor tip speed, flR, ft/sec 

v x 


increment in longitudinal velocity due to rotor inflow/wake, positive for 
ward, ft/sec 

v y 


increment in vertical velocity due to rotor inflow/wake, positive upward, 
ft/sec 

v z 


increment in sidewash velocity due to rotor inflow/wake, positive to the 
right when view is in direction of flight, ft/sec 

x,y,z 

(X,Y,Z) 

Cartesian coordinate system centered in rotor hub as shown in figure 1 , 
(convention of ref. 1 whereby x-z is the tip path plane and y upward, 
positive sense is as for the velocity increments, v x , v , and v z ) non- 
dimensionalized with rotor radius, /R 

r 

(GAMMA) 

blade circulation, ft 2 /sec 


(GAMSIN) 

factor for azimuthal variation of circulation as a function of sine of 
azimuth angle 


(GAMCOS) 

factor for azimuthal variation of circulation as a function of cosine of 
azimuth angle 


(GAMFACT) 

ratio of azimuthally affected circulation to assumed "average" circulation 

AA 

(DLAM) 

increment in rotor inflow ratio, v /Vj 

A n 

(DMU) 

increment in rotor advance ratio, v x /V t 


(VBAR) 

rotor advance ratio, V/V t 

A n 

(DNU) 

increment in sidewash velocity component, v z /V t 

p 

(RBAR) 

rotor blade station, r/R 

'a 


air density, slugs/ft 3 

V> 

(PSI) 

blade azimuth position, deg 

n 


rotor rotational speed, rad/s6c 
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CODE METHODOLOGY 


The development of the theory began with the consideration of a rotor vortex system 
which is formed mainly by shed vorticity leaving the rotor blades. The generation of this vor- 
ticity is a function of the radial variation of blade lift (i.e. "circulation"). As the horizontal 
flight speed increases the vortex column below the rotor deflects more to the rear. At suffi- 
ciently high speeds, the column becomes practically flat resulting in a relatively undistorted 
sheet, or ribbon of vorticity transported downstream at the free stream velocity. Vortices that 
constitute that ribbon (essentially attributable to the incremental changes in the radial variation 
of circulation) spring from all of the rotor blade sections. Thus, cycloidal patterns of vorticity 
of various shapes and dimensions are formed. In the analytical scheme and code these cycloidal 
vortices are represented by an equivalent rectilinear vortex system. References 2 and 3 offer a 
full description of the treatment of vorticity. 

In the computation of resultant induced velocities the contribution of lateral and lon- 
gitudinal vortices are computed separately and then summed. An important assumption in the 
theory development is that the circulation varies only with radius and that the contributions of 
azimuthal variations of vorticity are relatively small. However, a means for accounting for 
azimuthal variation is provided in the code and is described in Appendix A. 

The code is written in FORTRAN-77 and can be compiled on the Microsoft Version 5 
compiler. There are 500 lines of code which can be reduced to less than 350 lines if a user 
wishes to limit the calculation to a single field point velocity calculation (i.e. to be used as a 
subroutine, perhaps). The code is easily modified to the requirements of the user and requires 
no special hardware or software environment. A listing is given in Appendix B. 

The operation of the code begins with a selection of circulation (» lift) distribution on 
the rotor blade as a function of radius. The source code DOWN is compiled with one of three 
subroutines describing representative radial variations of circulation. The subroutines for several 
circulations as functions of blade radial station are identified as GA1, GA2, or GA3 and shown 
in Appendix C. GA1 is a linear distribution as shown on fig. 3.3 of reference 1. GA2 is a 
parabolic distribution described on page 56 of reference 1. GA3 is similar to equation 3.34 of 
reference 1. Within the code, the circulation is normalized so that the output velocities are nor- 
malized with tip speed. The "normalization" allows for any circulation pattern a user of the code 
may wish to try and is described in Appendix D. 

The organization of the code begins with the declaration of a parameter, NSEG, which 
is equated to the value of ten. NSEG is the number of blade segments for which incremental 
circulation will be calculated. NSEG can be increased for greater accuracy though there will be 
an associated increase in time required for calculation. Next, the user identifies the output file 
name and circulation subroutine called for in the compilation. The "normalization" factor for 
circulation is then calculated. Circulation (i.e. GAMMA) for the number of blade segments 
chosen is calculated and normalized. With that, the increment (or decrement) in circulation be- 
tween segments is determined. These will be used to calculate the velocity increments after the 
integration factors (K, L, M, N, which are described in reference 1) are completed. 

Before these integration factors are determined, in the main part of the program, the 
user chooses the location and type of velocity calculation. After the key input of circulation fac- 
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tors, advance ratio, and rotor thrust coefficient are given to the code, the dimensions required 
by the choice of calculation are input. These are the location of a plane, the initial, end, and 
incremental dimensions. That concludes the input and the code begins the primary loop ("DO 
50") for producing the integral factors K, L, M, and N. Before the leaving the loop the velocity 
increments effected by all blade segments are defined for a field point. These are then normal- 
ized i.e. converted to ratios of velocity increments divided by tip speed. 

USER INSTRUCTIONS 

The source code, DOWN (Appendix C), should be compiled with one of the chosen cir- 
culation subroutines (Appendix D) using a FORTRAN-77 compiler. The input to the executable 
file, then (from the keyboard in response to requests shown on the computer monitor), is : 

a. A name for the output file. 

b. Identification of circulation subroutine (GA1, GA2, GA3, or other). 

c. Choice of plane of calculation (calculation option 1, 2, 3, 4, or 5). 

1. Longitudinal- vertical (x-y) plane at lateral position z. 

2. Lateral- vertical (z-y) plane at longitudinal position x. 

3. Horizontal (x-z) plane at vertical height y. 

4. Radial variation at various azimuths (horizontal plane at height y). 

5. Azimuthal variation at various radii (horizontal plane at height y). 

d. Factors GAMSIN and GAMCOS accounting for azimuthal variation of circulation 

(described in Appendix A) 

e. Advance ratio, p. 

f. Rotor thrust coefficient, C T . 

Depending, then, on the choice (c) of plane of calculation (see figure 1) desired, input 
will be requested as follows: 

If 1 Lateral location of vertical plane, Z. 

Initial height, YINIT. 

Final height, YFINL. 

Incremental height, DELTAY. 

Initial longitudinal distance, XINIT. 

Final longitidinal distance, XFINL. 

Incremental longitudinal distance, DELTA X. 

If 2 Longitudinal location of vertical plane, Y. 

Initial lateral distance, ZINIT. 

Final lateral distance, ZFINL. 

Incremental lateral distance, DELTAZ. 

Initial height, YINIT. 

Final height, YFINL. 

Incremental height, DELTAY. 
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If 3 Height of horizontal plane, Y. 

Initial longitudinal distance, XINIT. 

Final longitudinal distance, XFINL. 

Incremental longitudinal distance, DELTAX. 

Initial lateral distance, ZINIT. 

Final lateral distance, ZFINL. 

Incremental lateral distance, DELTAZ. 

If 4 Height of horizontal plane, Y. 

Initial radial distance, RINIT. 

Final radial distance, RFINL. 

Incremental radial distance, DELTAR. 

Initial azimuth, AZINIT. 

Final azimuth, AZFINL. 

Incremental azimuth, DELTAP. 

If 5 Height of horizontal plane, Y. 

Initial azimuth, AZINIT. 

Final azimuth, AZFINL. 

Incremental azimuth, DELTAP. 

Initial radial distance, RINIT. 

Final radial distance, RFINL. 

Incremental radial distance, DELTAR. 

If "choice" c is 1, 2, or 3, then the coordinates X, Y, and Z of the field point are listed 
along with the incremental nondimensional velocities AA, An, and A rj and the corresponding 
radial position, p, and azimuthal coordinate, i>. The factor, GAMFACT, described in Appendix 
A, is listed as well. 

If "choice" c is 4 or 5, then the radial coordinate, p, and azimuthal coordinate, tp, are 
listed along with incremental nondimensional velocities AA, An, and At] and the factor, GAM- 
FACT. 


SAMPLE INPUT/OUTPUT 


Input 

Output file name, DwnDemo 

Circulation subroutine name used in compilation, GA3 

Calculation option, 5 

GAMSIN, 1.5 

GAMCOS, 1.12 

Advance ratio, VBAR, 0.149 

Rotor thrust coefficient, CT, 0.00630 
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Height, Y, 

0.077 

Initial azimuth, AZINIT, 

0 . 

Final azimuth, AZFINL, 

90. 

Increment in azimuth, DELTAP 

30. 

Initial radius, RINIT, 

0.2 

Final radius, RFINL, 

1.2 

Increment in radius, DELTAR 

0.2 


Output 

Output File, FNAME : DwnDemo 
Circulation program: GA3 

Radial variation at various Azimuths 

at Height, Y: .077; VBAR: .149; CT: .00630; GAMSIN: 1.500; GAMCOS: 1.120 


>1 

r/R 

DLAM 

DMU 

DNU 

GAMFACT 

0. 

.20 

-.00025 

.01301 

.02746 

1.0000 

0. 

.40 

-.00197 

.01757 

.03128 

1.0000 

0. 

.60 

-.00784 

.01987 

.03249 

1.0000 

0. 

.80 

-.01658 

.01896 

.03256 

1.0000 

0. 

1.00 

-.02665 

.00966 

.03121 

1.0000 

0. 

1.20 

-.02191 

.00168 

.03003 

1.0000 

30. 

.20 

-.01812 

.01172 

.02903 

.9007 

30. 

.40 

-.03214 

.01582 

.02172 

.9007 

30. 

.60 

-.04297 

.01789 

.01444 

.9007 

30. 

.80 

-.05509 

.01708 

.00454 

.9007 

30. 

1.00 

-.06681 

.00870 

-.01249 

.9007 

30. 

1.20 

-.05932 

.00151 

-.03101 

.9007 

60. 

.20 

-.02333 

.01098 

.01562 

.8437 

60. 

.40 

-.03120 

.01482 

.00555 

.8437 

60. 

.60 

-.03757 

.01676 

-.00673 

.8437 

60. 

.80 

-.04084 

.01600 

-.02364 

.8437 

60. 

1.00 

-.03639 

.00815 

-.05727 

.8437 

60. 

1.20 

.04101 

.00142 

-.04796 

.8437 

90. 

.20 

-.02165 

.01075 

.00815 

.8262 

90. 

.40 

-.02520 

.01452 

-.00254 

.8262 

90. 

.60 

-.02417 

.01642 

-.01465 

.8262 

90. 

.80 

-.01640 

.01565 

-.02870 

.8262 

90. 

1.00 

.01875 

.00776 

-.03536 

.8262 

90. 

1.20 

.01472 

.00138 

-.00496 . 

.8262 

120. 

.20 

-.01896 

.01098 

.00426 

.8437 

120. 

.40 

-.01900 

.01482 

-.00490 

.8437 

120. 

.60 

-.01400 

.01676 

-.01290 

.8437 
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330. 1.00 -.04961 .01086 .01241 1.1242 

330. 1.20 -.04605 .00189 .02620 1.1242 

DESCRIPTION of OPERATING ENVIRONMENT 

The source code "DOWN" was compiled with the Microsoft FORTRAN Version 5.0 on a 
Hewlett Packard VECTRA QS/20 (IBM compatible) computer. The Microsoft FORTRAN 
compiler operates on any IBM or IBM-compatible computer running MS-DOS Version 3.0 or 
later version. The operating system for the HP VECTRA was MS-DOS Version 4.0. The com- 
piled program, DOWN, can then be run under DOS Version 2.1 or later. The VECTRA had a 
math coprocessor chip and with its clock rate of 20 mhz a velocity calculation at a point was 
approximately 4 seconds. 
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APPENDIX A 


The azimuth variables, GAMSIN and GAMCOS, offer a means to account for variability 
of circulation with azimuth. According to reference 1 

T(p,n,rP) = T(p) + AIX/1,0) 

Though the theory is developed with AT(p.,i>) « 0.0, azimuth dependency may be introduced as 
on page 95 of reference 1: 

Ar(/x,0) » -1.5 p sin ip + 1.12 /i 2 (1 - cos 2 $) 


In the code then: 

Ar(/i,V>) » -(GAMSINXVBAR) sin(PSI) + (GAMCOS)(VBAR)(VBAR)(l - cos(2*PSI)) 
where GAMSIN = 0.0 or 1.5 

and GAMCOS = 0.0 or 1.12 

In the output of the code an accounting for azimuth dependency is given as GAMFACT where: 
GAMFACT = ( r(p) + Ar(/i,V0) / r (p) 
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APPENDIX B 


PROGRAM DOWN 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

C CALCULATION OF ROTOR DOWNWASH FIELD USING FLAT-WAKE THEORY. 

C PROGRAM ORIGINALLY WRITTEN BY M. HAGLUND, ’91, PRINCETON UNIV. 

C MODIFIED BY R. MCKILLIP 7/90 FOR IBM-PC USAGE. 

C (ADDITIONAL MODIFICATIONS BY J. WILSON 10/90 — 2/91) 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
NSEG IS NUMBER OF ROTOR SEGMENTS 

IMPLICIT INTEGER (I-N) 

IMPLICIT REAL (A-H,0-Z) 

REAL J,K,KINF,L,LAMBDA,LINF,M,MINF,MU,N,NINF,NU 
PARAMETER (NSEG=10) 

C 

REAL GAMMA(NSEG),DGAMMA(NSEG),DVYI(NSEG),DVZI(NSEG),DVXI(NSEG) 
CHARACTER* 10 FNAME 
CHARACTER* 10 CIRC 
C 

DATA PI/3.141592654/ 

C 

100 FORMAT( 1 X,79A) 

WRITE(* 100) ’ ****************************************’ 

WRITER, 100) ’ * FLAT WAKE DOWNWASH CALCULATION *’ 

WRITER 100) ’ ***************************************** 

D2R = PI/ 180. 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

C OPEN FILES FOR OUTPUT 

C UNIT 1 IS DATA; PARAMETERS ARE PRINTED ON TITLE LINE. 

C 

WRITE(*,100) ’ Enter filename for output data:’ 

READ(*,’(A)’) FNAME 
OPEN( 1 ,FILE=FNAME) 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

WRITE(1,8) FNAME 

8 FORMAT(lX,’ Output File, FNAME :’,10A) 

WRITE(*,100) ’ Enter Circulation Program used (i.e. GA1, etc.) ’ 

READ(*,’(A)’) CIRC 

WRITE(1,9) CIRC 

9 FORMAT(lX,’ Circulation program ’,10A) 
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c 

C CIRCULATION GAMMA IS COMPUTED USING A FUNCTION OF THE FORM: 

C GAMMA = GA(RHO), WHICH IS COMPILED SEPARATELY AND LINKED TO 

C THIS CODE, AS APPROPRIATE FOR THE TARGET COMPUTER. 

C 

C THIS SECTION COMPUTES BOTH THE CIRCULATION GAMMA(I) AND THE 

C INCREMENTAL CIRCULATION DGAMMA(I). 

C 

C FIRST, CALCULATE THE NORMALIZATION FACTOR 

C 

G = 0.0 
NSTEPS = 200 
BG = 0.0 
ED = 1.0 

DRHO = (ED-BG)/FLOAT(NSTEPS) 

RHO = BG 

RH02 = RHO + DRHO 
GAMA = GA(RHO) 

4 CONTINUE 

IF (RH02.LT.ED) THEN 
GAMMAD = GA(RH02) 

G = G + (GAMA*RHO) + (GAMMAD*RH02) 

GAMA = GAMMAD 
C 

RHO = RH02 
RH02 = RH02 + DRHO 
GOTO 4 
ENDIF 

G = G * DRHO 
C 

C NOW CALCULATE THE CIRCULATION AND NORMALIZE BY G 

C 

DO 5 I=1,NSEG 

RHO = ,5*(2*I-l)/FLOAT(NSEG) 

GAMMA(I) = GA(RHO)/G 

5 CONTINUE 
GAMMA(l) = 0.0 

C THIS LINE HAS BEEN REPLACED BY THE DO LOOP THAT FOLLOWS. 

C GAMMA(Il) = 0.0 

C 

C ...AND CONVERT NORMALIZED CIRCULATION TO A STAIR-STEP 

C FUNCTION FOR THE CALCULATION. 

C 
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DO 6 1 = 1, NSEG-1 

DGAMMA(I) = GAMMA(I) - GAMMA(I+1) 

6 CONTINUE 
DGAMMA(NSEG) = GAMMA(NSEG) 

CCCCCCCCC 1 -CCCCCCCC2-CCCCCCCC3-CCCCCCCC4-CCCCCCCC5-CCCCCCCC6- CC 
C 

C ASK FOR LOCATION AND TYPE OF VELOCITY CALCULATION 

C 

7 CONTINUE 

WRITE(*,100) ’ Enter Calculation Option:’ 

WRITE(*,100) ’ 1) Longitudinal/Vertical plane at Lateral dist.’ 

WRITE(*,100) ’ 2) Lateral/ Vertical plane at Longitudinal dist.’ 

WRITE(*,100) ’ 3) Horizontal plane at Vertical Height’ 

WRITE(*,100) ’ 4) Radial variation at various Azimuths’ 

WRITE(*,100) ’ 5) Azimuthal variation at various Radii’ 

WRITER, 100) ’ ?: ’ 

C 

READ(V) IDIR 

IF( (IDIR.LT. 1 ).OR.(IDIR.GT.5) ) GO TO 7 
WRITE(*,100) ’ Enter factors for circulation varying with az.’ 

READ(*,*) GAMSIN,GAMCOS 
C 

WRITE(*,100) ’ Enter ADVANCE RATIO: ’ 

READ(V) VBAR 
C 

WRITE(MOO) ’ Enter ROTOR THRUST COEFFICIENT: ’ 

READ(V) CT 
C 

C INCREMENT ALONG X- DIRECTION (and Y) for long./vert. plane 

C 

IF( IDIR.EQ.l ) THEN 

WRITE(*,100) ’ Enter lateral location of Vertical plane:’ 

READ(V) Z 

WRITE(MOO) ’ Enter initial, final, and delta Height’ 

READ(*,*) YINIT,YFINL,DELTAY 

WRITE(*,100) ’ Enter initial, final, and delta Long, dist.’ 

READ(*,*) XINIT,XFINL,DELTAX 

WRITE(1,100) ’ Longitudinal/Vertical plane at Lateral dist.’ 

WRITE( 1,10) Z, VB AR,CT,G AMSIN,G AMCOS 
10 FORMAT(lX,’ at, Z:’,F6.3,’; VBAR:’,F5.3,’; CT:’, 

> F7.5,’; GAMSIN:’,F5.3,’; GAMCOS:’,F5.3) 

WRITE(1,100) ’ X Y Z DLAM DMU DNU 

> GAMMA R PHI’ 
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WRITE(1,100) ’ 


X = XINIT - DELTAX 
Y = YINIT 

INCREMENT ALONG Y- DIRECTION (and Z) for lateral/vertical plane 

ELSE IF ( IDIR.EQ.2 ) THEN 

WRITE(*,100) ’ Enter Long, location of Vert, plane’ 

READ(V) X 

WRITE(*,100) ’ Enter init., final, and delta Lateral pos.’ 

READ(V) ZINIT,ZFINL,DELTAZ 
WRITE(MOO) ’ Enter init., final, and delta Height ’ 

WRITE(1,100) ’ Lateral/ Vertical plane at Longitudinal dist.’ 
READ(*,*) YINIT, YFINL,DELT AY 
WRITE( 1,11) X, VB AR,CT,G AMSIN,GAMCOS 
F0RMAT(1X,’ at, X:’,F6.3,’; VBAR:’,F5.3,’; CT:\ 

F7.5,’; GAMSIN:’,F5.3,’; GAMCOS:’,F5.3) 

WRITE(1,100) ’ X Y Z DLAM DMU DNU 
GAMMA R PHI’ 

WRITE(1,100) ’ - 


Z = ZINIT - DELTAZ 
Y = YINIT 

INCREMENT ALONG Z-DIRECTION (and X) for horizontal plane 

ELSE IF ( IDIR.EQ.3 ) THEN 

WRITE(*,100) ’ Enter height of Horizontal plane’ 

READ(*,*) Y 

WRITE(*,100) ’ Enter init., final, and delta Horiz. distance’ 

REAEKV) XINIT, XFINL,DELTAX 

WRITE(*,100) ’ Enter init., final, and delta Lateral distance’ 

READ(*,*) ZINIT, ZFINL, DELTAZ 

WRITE(1,100) ’ Horizontal plane at Vertical Height’ 

WRITE( 1,12) Y, VB AR,CT,G AMSIN,G AMCOS 
F0RMAT(1X,’ at, Y:’,F6.3,’; VBAR:’,F5.3,’; CT:’, 

F7.5,’; GAMSIN:’,F5.3,’; GAMCOS:’,F5.3) 

WRITE(1,100) ’ X Y Z DLAM DMU DNU 
GAMMA R PHI’ 

WRITE(1,100) ’ — 


X - XINIT - DELTAX 



Z = ZINIT 
C 

C INCREMENT ALONG RADIAL - for various azimuths for a horiz. plane 

C 

ELSE IF ( IDIR.EQ.4 ) THEN 

WRITE(*,100) ’ Enter Height of horizontal plane’ 

READ(V) Y 

WRITE(*,100) ’ for Radial variation at various Azimuths’ 
WRITE(*,100) ’ Enter initial, final, and delta Radius’ 

READ(V) RINIT,RFINL,DELTAR 

WRITE(*,100) ’ Enter initial, final, and delta Azimuth’ 

READ(V) AZINIT,AZFINL,DELTAP 
WRITE(1,100) ’ Radial variation at various Azimuths’ 

WRITE( 1,13) Y,VB AR,CT,G AMSIN,GAMCOS 
13 F0RMAT(1X,’ at Height, Y:’,F6.3,’; VBAR:’,F5.3,’; CT:’, 

> F7.5,’; GAMSIN:’,F5.3,’; GAMCOS:’,F5.3) 

WRITE(1,100) ’ AZ R DLAM DMU DNU 

> GAMMA’ 

WRITE(1,100) ’ 


RBAR = RINIT - DELTAR 
PHI = AZINIT 

C 

C INCREMENT ALONG ANNULUS for various radii for a horiz. plane 

C 

ELSE IF (IDIR.EQ.5) THEN 

WRITE(*,100) ’ Enter height of horizontal plane’ 

C 

READ(*,*) Y 

WRITE(*,100) ’ for Azimuthal variation at various Radii’ 

WRITE(*,100) ’ Enter initial, final, and delta Azimuth’ 

READ(*,*) AZINIT, AZFINL,DELTAP 
WRITER, 100) ’ Enter initial, final, and delta Radii’ 

READ(*,*) RINIT,RFINL, DELTAR 

WRITE( 1,100) ’ Azimuthal variation at various Radii’ • 

WRITE(1,14) Y,VBAR,CT,GAMSIN,GAMCOS 
14 F0RMAT(1X,’ at Height, Y:’,F6.3,’; VBAR:’,F5.3,’; CT:’, 

> F7.5,’; GAMSIN:’,F5.3,’; GAMCOS:’,F5.3) 

WRITE(1,100) ’ AZ R DLAM DMU DNU 

> GAMMA’ 

WRITE(1,100) ’ 


RBAR = RINIT 
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PHI = AZINIT - DELTAP 
AZFINLX = AZFINL - DELTAP 
ENDIF 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

LOOP POINT FOR INTEGRATION CALCULATIONS 

WRITE(*,’( 1 X,80 A)’) ’ CALCULATING. . 

15 CONTINUE 

IF ( (IDIR.EQ.1).AND.(X.GE.XFINL).AND.(Y.GE.YFINL) ) 

> GOTO 999 

IF ( (IDIR.EQ. 1 ). AND.(X.GT.XFINL) ) THEN 
X = XINIT 

Y = Y + DELTAY 
ELSE IF (IDIR.EQ. 1) THEN 

X = X + DELTAX 
ENDIF 
C 

IF ( (IDIR.EQ.2).AND.(Z.GE.ZFINL).AND.(Y.GE.YFINL) ) 

> GOTO 999 

IF ( (IDIR.EQ.2).AND.(Z.GT.ZFINL) ) THEN 
Z = ZINIT 

Y = Y + DELTAY 
ELSE IF (IDIR.EQ.2) THEN 

Z = Z + DELTAZ 
ENDIF 
C 

IF ( (IDIR.EQ.3).AND.(X.GE.XFINL).AND.(Z.GE.ZFINL) ) 

> GOTO 999 

IF ( (IDIR.EQ.3).AND.(X.GT.XFINL) ) THEN 
X = XINIT 
Z = Z + DELTAZ 
ELSE IF (IDIR.EQ.3) THEN 
X = X + DELTAX 
ENDIF 
C 

IF ( (IDIR.EQ.4).AND.(RBAR.GE.RFINL). AND. (PHI.GE. AZFINL) ) 

> GOTO 999 

IF ( (IDIR.EQ.4).AND.(RBAR.GT.RFINL) ) THEN 
RBAR = RINIT 
PHI = PHI + DELTAP 
ELSE IF (IDIR.EQ.4) THEN 
RBAR = RBAR + DELTAR 


15 


ENDIF 


IF ( (IDIR.EQ.5).AND.(PHI.GE.AZFINL).AND.(RBAR.GE.RFINL) ) 

GOTO 999 

IF ( (IDIR.EQ.5).AND.(PHI.GT.AZFINLX) ) THEN 
PHI = AZINIT 
RBAR = RBAR + DELTAR 
ELSE IF (IDIR.EQ.5) THEN 
PHI = PHI + DELTAP 
ENDIF 

IF (IDIR.LE.3) THEN 

RBAR = SQRT (X*X + Z*Z) 

IF (RBAR.LT.0.01) RBAR = .01 
SINAZ = Z/RBAR 
COSAZ = -X/RBAR 
AZ * ASIN(SINAZ) 

PHI = AZ/D2R 

IF ((SINAZ.GE.O.).AND.(COSAZ.GE.O.)) PHI = PHI 
IF ((SINAZ.GE.O.).AND.(COSAZ.LT.O.)) PHI = 180. - PHI 
IF ((SINAZ.LT.O.).AND.(COSAZ.LT.O.)) PHI = 180. - PHI 
IF ((SINAZ. LT.O.). AND. (COSAZ.GE.O.)) PHI = 360. + PHI 
COS2AZ = COS(2*PHI*D2R) 

GAMFACT = 1. - GAMSIN*VBAR*SINAZ + GAMCOS*VBAR*VBAR*(l. 
COS2AZ) 

ENDIF 

IF ( IDIR.GE.4 ) THEN 

X = -RBAR*COS(PHI*D2R) 

Z = RBAR*SIN(PHI*D2R) 

SINAZ = SIN(PHI*D2R) 

COS2AZ = COS(2*PHI*D2R) 

GAMFACT = 1. - GAMSIN*VBAR*SINAZ + GAMCOS*VBAR*VBAR*(l. 
COS2AZ) 

ENDIF 

VX = 0.0 
VY = 0.0 
VZ = 0.0 

DO 50 II = 1 ,NSEG 


RHO = II/FLOAT(NSEG) 


VBARSTAR = VBAR/RHO 

XI = X/RHO 

XINF - -20.0/RHO 

Y1 = Y/RHO 

YINF = Y1 

Z1 = Z/RHO 

ZINF = Z1 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

C FIRST INTEGRAL AROUND DISK FROM PI/2 TO 3PI/2, PRODUCING 

C FACTORS K,L,M AND N. BG AND ED ARE THE LIMITS OF INTEGRATION, 

C S IS INTEGRAND. 

C 

BG = PI/2. 

ED - 3.*PI/2. 

S = 0.01 

SINV2P = S / (2.0*PI) 

A = BG 
K =■ 0.0 
KINF = 0.0 
L = 0.0 

LINF = 0.0 
M = 0.0 
MINF = 0.0 
N = 0.0 

NINF = 0.0 
U = 0.0 

UINF = 0.0 
B1 = 0.0 
B1INF = 0.0 
B2 = 0.0 
C 

B2INF = 0.0 
W1A = 0.0 
W1AINF = 0.0 
W2A = 0.0 
W2AINF = 0.0 
W1B = 0.0 
W1BINF = 0.0 
W2B = 0.0 
W2BINF = 0.0 
W1 = 0.0 
W1INF = 0.0 
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non 


W2 = 0.0 
W2INF = 0.0 


C 

20 CONTINUE 

IF (A.LT.ED) THEN 
CS1 = COS(A) 

CS2 = COS(A+S) 

SN1 = SIN(A) 

SN2 = SIN(A+S) 

B1 = Y1*Y1 + (SN1-Z1)*(SN1-Z1) 

B1INF = YINF*YINF + (SN1-ZINF)*(SN1-ZINF) 

B2 = Y1*Y1 + (SN2-Z1)*(SN2-ZI) 

B2INF = YINF*YINF + (SN2-ZINF)*(SN2-ZINF) 

W1A = (CS1 -X 1 )/SQRT((CS 1 -X 1 )*(CS1 -X 1 )+B 1 ) 

W1AINF = (CS 1 -XINF)/SQRT((CS1 -XINF)*(CS 1 -XINF)+B 1 INF) 

W1B = (CS 1 +X 1 )/SQRT((CS 1 +X 1 )*(CS 1 +X 1 )+B 1 ) 

W1BINF = (CS 1 +XINF)/SQRT((CS 1 +XINF)*(CS1 +XINF)+B 1 INF) 

W1 = W1A - W1B 
W1INF = W1AINF - W1BINF 
W2A = (CS2-X 1 )/SQRT((CS2-X 1 )*(CS2-X 1 )+B2) 

W2AINF = (CS2-XINF)/SQRT((CS2-XINF)*(CS2-XINF)+B2INF) 

W2B = (CS2+X 1 )/SQRT((CS2+X 1 )*(CS2+X 1 )+B2) 

W2BINF = (CS2+XINF)/SQRT((CS2+XINF)*(CS2+XINF)+B2INF) 

W2 = W2A - W2B 
W2INF = W2AINF - W2BINF 
K. = K. + (Wl/Bl) + (W2/B2) 

KINF = K.INF + (W1INF/B1INF) + (W2INF/B2INF) 

L = L + (SN1*W1/B1) + (SN2*W2/B2) 

LINF = LINF + (SN 1 * W 1 INF/B 1 INF) + (SN2*W2INF/B2INF) 

U = U + (SN1*SN1*W1/B1) + (SN2*SN2*W2/B2) 

UINF = UINF + (SN 1 *SN 1*W1 INF/B 1 INF) + (SN2*SN2*W2INF/B2INF) 
A = A + S 
GO TO 20 
ENDIF 

SWAP ORDER OF CALCULATION TO AVOID 0.0/0.0 SINGULARITY 

N = (U-Z1*L)*SINV2P 
NINF = (UINF-Z1*LINF)*SINV2P 
M = (L-Z1*K)*SINV2P 
MINF = (LINF-Z1*KINF)*SINV2P 
C 

K = K*SINV2P*Y1 
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KINF = KINF*SINV2P*YINF 

L = L*SINV2P*Y1 

LINF = LINF*SINV2P*YINF 

U = U*SINV2P 

UINF = UINF*SINV2P 

C (OLD METHOD USED UPDATED K,L,U VALUES) 

C N = U-(Z1*L/Y1) 

C NINF = UINF-(ZINF*LINF/YINF) 

C M = (L-Z1*K)/Y1 

C MINF = (LINF-ZINF*KINF)/YINF 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c 

C SECOND INTEGRAL OVER LIMITS -1 TO 1 PRODUCES FACTORS H AND J. 

C 

H = 0.0 
HI = 0.0 
H1A = 0.0 
H2B = 0.0 
H2 = 0.0 
H2A = 0.0 
H2B = 0.0 
YX = 0.0 
YX1 1 = 0.0 
YX21 = 0.0 
YX31 = 0.0 
YX12 = 0.0 
YX22 = 0.0 
YX32 = 0.0 
BG = -1.0 
ED = 1.0 
A = BG 
C 

YX1 1 = Y1/(Y1*Y1 + (A-Z1)**2) 

YX21A = XI + SQRT(1-A*A) 

YX21B = SQRT((X1 + SQRT(1-A*A))**2 + Y1*Y1 + (A-Z1)**2) 

YX21 = YX21 A/YX21B 
YX31A = XI - SQRT(1-A*A) 

YX31B = SQRT((X1 - SQRT(1-A*A))**2 + Y1*Y1 + (A-Z1)**2) 

YX31 = YX31A/YX31B 
YX01 = YX1 1*(YX21 - YX31) 

H1A = 1./YX21B 
H1B = 1/YX31B 
HI = HI A - H1B 
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CONTINUE 


30 
C 

IF (A.LE.(ED-S)) THEN 
A = A + S 

YX12 = Y1/(Y1*Y1 + (A - Zl)**2) 

YX22A = XI + SQRT(1-A*A) 

YX22B = SQRT((X1 + SQRT(1-A*A))**2 + YI*Y1 + (A-Z1)**2) 

YX22 = YX22A/YX22B 
YX32A = XI - SQRT(1-A*A) 

YX32B = SQRT((X1 - SQRT(1-A*A))**2 + Y1*Y1 + (A-Z1)**2) 

YX32 = YX32A/YX32B 
YX02 = YX12*(YX22 - YX32) 

YX = YX + YX01 + YX02 
YX01 = YX02 
H2A = 1./YX22B 
H2B = 1./YX32B 
H2 = H2A - H2B 
H = H + HI + H2 
HI = H2 
GOTO 30 
ENDIF 

H = H * (-1.) * 0.5 * SINV2P 
J = YX * (-1.) * SINV2P * 0.5 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

C INTEGRAL FACTORS COMPLETE. NOW BUILD INFLUENCE COEFFS, AND 
C MULTIPLY BY CIRCULATION DISTRIBUTION TO GET VELOCITIES. 

C 

DVYI(II) = H - 0.5*(VBARSTAR * (MINF + M) + (NINF + N)) 

DVZI(II) = -0.5 * (VBARSTAR * (KINF + K) + (LINF + L)) 

DVXI(II) = J 
C 

VY = VY + DVYI(II) * DGAMMA(II) * GAMFACT 
VZ = VZ + DVZI(II) * DGAMMA(II) * GAMFACT 
VX = VX + DVXI(II) * DGAMMA(II) * GAMFACT 
50 CONTINUE 

C 

C NORMALIZE THE VELOCITY COMPONENTS AND PRINT OUTPUT 

C 

VO = CT / (2.0*VBAR) 

LAMBDA = VY * VO 
MU = -1.0 * VX * VO 
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NU = VZ * VO 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccc cccccccc 
c 

C UPDATE DESIRED EVALUATION POSITION, AND BRANCH BACK 

C 

IF (IDIR.LE.3) THEN 

WRITE(1,56) X, Y, Z, LAMBDA, MU, NU, GAMFACT,RBAR,PHI 
WRITE(*,56) X, Y, Z, LAMBDA, MU, NU, GAMFACT,RBAR,PHI 

56 FORMAT(1X,3(F6.2,1X),3(F10.5,1X),F8.4,F6.2,F8.0) 

ELSE 

WRITE(1,57) PHI, RBAR, LAMBDA, MU, NU, GAMFACT 
WRITE(*,57) PHI, RBAR, LAMBDA, MU, NU, GAMFACT 

57 FORMAT(1X,F8.0,2X,F6.2,2X,3(F10.5,2X),F9.4) 

ENDIF 

C 

GO TO 15 
C 

C END OF ITERATIONS; CLOSE DATA FILE. 

C 

999 CONTINUE 

CLOSE(l) 

STOP 

END 
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APPENDIX C 


CIRCULATION SUBROUTINES 


C 

C GA1 FUNCTION FOR THE CIRCULATION DISTRIBUTION 
C Linear distribution (fig. 3.3 of ref. 1) 

C 

REAL FUNCTION GA(A) 

REAL A 
C 

GA = 10.0 * (0.0589 + 1.3783* A) 

C 

RETURN 

END 


C 

C GA2 FUNCTION FOR THE CIRCULATION DISTRIBUTION 

C Parabolic distribution (page 56 of ref. 1) 

C 

REAL FUNCTION GA(A) 

REAL A 
C 

GA = 10.0 * (A* A - A**3) 

C 

RETURN 

END 


C 

C GA3 FUNCTION FOR THE CIRCULATION DISTRIBUTION 

C EQ. 3.34 (of ref. 1) for blades with -10 deg twist 

C 

REAL FUNCTION GA(A) 

REAL A 
C 

GA = ((1.0 - 0.685*A)*A*(A*A + 0.011))/(A*A + 0.0055) 

C 

RETURN 

END 
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APPENDIX D 


In the code the circulation is normalized, so that any pattern can be coupled to the 


code, as follows: 


G(normal) s 2 1 T(p) p dp 

0 


2 / 

J 0 


1 


m 

r(normal) = = — -- 


r(p) 

G 


2 / r(p) p dp 

Jo 

This results in velocity components normalized by the momentum value of the downwash inde- 
pendent of the choice of T(p), where : 

v y A i 

Vy = -- = 

V V 

0 v o 


and where 


v z 


V V 

o v o 


v x -I*i 

v„ = — = — 


V V 

0 o 


v 0 C T 


v o = 


2p 


and then as final output 
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